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Summary 


In  present  investigation,  analyses  for  the  damage  evolution  behavior  of  particulate 
composite  materials  by  using  the  finite  element  method  (FEM)  and  the  s-version  finite 
element  method  (s-FEM)  were  carried  out.  The  analyses  were  carried  out  in  particular 
interest  in  the  phenomenon  of  crack  propagation. 

Prior  to  crack  propagation,  material  damage  develops  in  the  material.  The  material 
damage  may  be  in  the  forms  of  microviod  and/or  microcracks  in  the  binder  (matrix)  and 
in  the  form  of  binder  (matrix)/particle  separation  that  is  known  to  be  dewetting.  In  a 
macroscopic  sense,  the  reinforcing  particles  distribute  evenly  in  matrix.  However,  at 
microscopic  level,  the  density  of  the  distributed  particles  varies.  This  means  that  the 
stiffness  and  strength  of  the  material  also  have  some  spatial  variations.  Material 
damages  initiate  at  the  weak  material  locations  and  then  propagate  the  surroundings. 
When  cracks  are  present  in  the  material,  the  cracks  interact  with  the  surroundings  and 
the  material 

To  simulate  such  scenarios,  we  adopted  two  kinds  of  damage  constitutive  models.  One 
is  isotropic  damage  model  and  the  other  is  “separate  dilatational/deviatoric  damage 
constitutive  model”  in  which  the  contributions  of  hydrostatic  and  of  deviatoric  stresses 
are  accounted  for  independently.  A  parameter  in  the  separate  dilatational/deviatoric 
damage  model  can  characterize  which,  hydrostatic  or  deviatoric  stress  component,  has 
dominant  influence  to  the  damage  behavior  of  the  material. 

A  series  of  analyses  on  uncracked  and  cracked  specimen  with  statistically  varying 
material  stiffness  at  a  microscopic  level  were  carried  out.  The  results  revealed  that  the 
damage  behavior  is  highly  influenced  by  the  damage  mode. 


3 


1.  Damage  constitutive  law 


In  present  investigation,  we  adopt  the  continuum  damage  constitutive  law  of  Simo  and 
Ju  [1,  2],  The  damage  model  was  extended  so  that  it  can  model  the  isotropic,  deviatoric, 
dilatational  and  deviatoric/dilatational  combined  damages  was  adopted.  The  constitutive 
model  was  implemented  in  the  in-house  s-FEM  computer  program.  The  s-FEM 
computer  program  of  Okada  et  al.  [3]  can  model  particulate  composite  materials  by 
superposing  the  finite  element  models  for  the  particles  on  that  for  the  global  structure  or 
for  the  region  of  unit  cell.  A  finite  element  model  for  a  particle  and  its  surrounding 
material  region  needs  to  be  built  only  once.  Such  finite  element  model  is  called  the 
“local  model”.  Then  it  is  repeatedly  superposed  on  the  one  that  represents  the  global 
structure,  that  is  called  the  “global  model”. 

In  present  investigation,  an  ordinary  finite  element  method  is  mostly  used.  The  s-FEM 
computer  program  can  also  be  used  to  carry  out  ordinary  finite  element  analyses.  To  do 
so,  we  simply  do  not  superpose  any  local  finite  element  models  on  the  global  one. 

In  this  chapter,  we  first  explore  the  damage  constitutive  model  that  is  adopted  in  this 
research. 

1.1  Isotropic  damage 

First,  we  describe  the  isotropic  damage  theory  by  following  Simo  and  Ju  [1],  In  Simo 
and  Ju  [1],  the  effective  stress  concept  and  the  hypothesis  of  strain  equivalence  are 
described.  When  they  are  applied  to  the  case  of  elastic  damage,  the  elastic  potential 

energy  of  damaged  material,  in  terms  of  the  strains  s and  the  damage  parameter  d 
is  written  to  be: 


y/(eu,d)  =  {\-d}f/0{£U) 


(1) 


where  ¥°(£ke)  is  the  elastic-potential  function  for  virgin  material  that  is  written  to  be: 


¥ 


1%)  =  lCi- 


ijkl  £ij  skl 


(2) 


Just  like  equivalent  stress  concept  in  the  theory  of  plasticity,  a  scalar  parameter  r  that 
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measures  the  magnitude  of  stress/deformation  is  introduced. 


r  =A/2^°(fW)  (3) 

Criteria  for  the  damage  evolution  are  written  as  follows.  First,  the  effective  stress  needs 
to  reach  the  critical  value. 

g{* (£k( )  -  r(d ))  =  f(eke  )-r{d)  =  0  (4) 

where  r(d )  is  the  function  of  the  damage  parameter  d  .  The  relationship  between  r(d) 
and  d  needs  to  be  determined  based  an  experimental  data.  The  material  damage 
progresses  when  the  damage  parameter  increases.  Therefore,  when  the  damage  is 
ongoing,  the  time  derivative  of  the  damage  parameter  d  must  be  positive. 

d  > 0  (5) 

d  is  determined  through  the  evolution  law,  as: 

d  =  rH(f,d)=  i{sy  )h(t,  d )  (6) 

The  constant  H(f,d)  characterizes  the  progress  of  damage  with  respect  to  the  effective 
stress  like  term.  Equations  (4)  and  (5)  must  be  satisfied  when  the  progressive  damage 
takes  place.  It  is  noted  that  the  constant  H(f,d )  is  always  positive.  Therefore,  the 
equivalent  stress-like  term  must  increase  to  satisfy  equation  (6),  while  the  material 
damage  is  in  progress. 

Since  the  elastic  potential  energy  is  assumed  in  the  form  of  equation  (1),  the  stresses 
a ij  are  written  in  the  following  form. 


aij  d)Cijkl£ld 


(7) 


The  rates  of  stresses  can  be  expressed  by  differentiating  both  the  sides  of  equation  (7), 
as: 
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^  ij  (l  dY-ijkfZ  kt:  dCijfcgSfcg  (l  d]C,jkfSkf  d(7jj 


(8) 


where  a ij(=Cijk(£k()  are  the  stresses  when  the  virgin  material  is  assumed  for  the  same 
strains.  From  equation  (6),  the  rate  d  of  damage  parameter  d  can  be  derived  to  be: 


d  =  tH(r,d)=  H{t  ,d)(y^r  e  k( 

T 


(7ke£ke 


(9) 


H(z,d )  characterizes  the  rate  of  damage  parameter  d  with  respect  to  the  rate  of  the 
effective  stress-like  parameter  t . 

Substituting  equation  (9)  in  equation  (8),  we  arrive  at: 


°ij  ~  (l _  d)Cijke£k(  -  dcr°j  -  (l  -d]CijU£ke - ~  ~ <rij<Tu*kt 


~  Dijk(<£kf 


(10) 


where 

Df,=(l (11) 
are  the  tangent  moduli  for  the  isotropic  damage  material  and  have  the  major  and 
minor  symmetries,  i.e.  ofkt  =  ,  D$t  =  D$k  and  D,™  =  D% . 


1.2  Separate  Isotropic/deviatoric  damage  model 

The  damage  evolution  of  some  class  of  materials  are  more  sensitive  to  hydrostatic 
pressure  stress  than  shear  stresses.  A  typical  scenario  is  in  a  material  containing 
microvoids.  Microvoids  grow  under  applied  positive-hydrostatic  pressure  stress.  But 
they  do  not  grow  under  negative-hydrostatic  pressure  stress.  The  growths  of  the  voids 
are  assumed  to  be  less  sensitive  to  the  shear  (deviatoric)  stresses  than 
positive-hydrostatic  stress. 
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One  way  to  model  such  material  is  to  separate  the  contributions  of  hydrostatic  and 
deviatoric  stresses  to  the  damage  growth.  To  do  so  with  a  very  simple  model,  we 
separate  the  damage  parameter  d  into  the  dilatational  (volumetric)  and  deviatoric  parts. 
Hence,  the  elastic  potential  energy  function  [equation  (1)]  is  modified  and  is  written  to 
be: 


y/{ski ,dv,dD)-(l-dv  (ekk )  +  (l  -  dD Vd 


(12) 


where  ekk  and  e-  are  the  volumetric  and  deviatoric  strains.  The  functions  Yv  and 
Y°d  are  defined  to  be: 

Yv=^Ki£kk)2  and  Yv  =  ^44  (13) 

where  K  and  //  are  the  bulk  and  shear  moduli.  It  is  noted  here  that  under  a  constraint 
condition  dv  =  dD ,  the  separate  dilatational/deviatoric  damage  model  is  the  same  as  the 
isotropic  damage  model.  We  define  two  kinds  of  effective  stress-like  terms,  as: 

=  ^Yv  (£kk  )  =  yjK(£kk  )2  and  =  tJ^-Yo  (£ij  )  =  (14) 

When  the  damages  are  in  progress,  the  following  conditions  must  be  satisfied 
Ty  =  rv  (dv )  and  akk  >  0  for  the  dilatational  damage 
fD  =  rv(dD)  for  the  deviatoric  damage 
The  evolution  equations  for  the  damage  parameters  are  written,  as: 

Ks 

dy  =  TyHy{fy  ,dy)=  Hy(fy  ,dvy/K£kk  =  Hy  (fy  ,  dy)—j===  £kk  (17) 

aHJ2 


(15) 

(16) 
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(18) 


d  ft  ~tdH  d(tD’^d)~  H  d{td^d) 


2{l£'kf 


pq£  pq 


r£ld 


The  constants  Hv(fv,dv )  and  HD(fD,dD )  govern  the  evolution  laws  for  the  damage 
parameters.  The  constants  are  the  functions  of  the  effective  stress-like  terms  fv  and 
fD  and  of  the  damage  parameters  dv  and  dD .  From  equation  (12),  we  can  write  the 
stresses,  as: 


<7;; 


dv{£kl’dv  >dD) 


?)£:, 


-  (l  ~dv  )djj  Kskk  +  2(l  -  dD  )//£ 


(19) 


By  differentiating  both  the  sides  of  equation  (19)  and  making  use  of  equations  (17)  and 
(18),  we  can  write  the  rate  form  constitutive  equation,  as: 


-  (l  ~  dy]SijKskk  +  2(l  -  d  D  -dvKskk  -2dDne\j 

=  (l-dv ]S{j Kakk  +  2(l -dD )//£■,)  - S,j  — (Kskk )2 su  —  ^~D ’  - ^ // 2 s-j s'u £ke 

TV  ZD 

-  r>v~Dp 
~Uijtd  £kl 


and, 


(20) 


D 


lu  =  (1  ~dv  )KSySu  +  (1  -  dD  )ju(sikSjt  + 

TV  tD 


djk  d;t ) 
M2£ij£'u 


(21) 


Dyke°  are  the  tangent  moduli  that  relate  the  rate  of  stresses  to  those  of  strains.  It  is  noted 


that  D;jud  have  the  major  and  minor  symmetries,  i.e.  Dyke°  =  DklijD ,  Dyke°  =  D\Jfk  and 


dv~d  =  dv-d 
u ijki  jikl 


1.3  Determining  damage  evolution  law  for  the  isotropic  damage  material 

The  damage  evolution  law  of  equation  (9)  has  a  scalar  function  H(f,d)  of  f  and  d  . 
H(f,d)  characterizes  the  damage  evolution  behavior  and  should  be  derived  from  a  set 
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of  experimental  data.  In  this  section,  how  the  scalar  function  H(r,d)  is  determined  is 
described. 

First,  we  assume  that  we  have  a  stress-strain  curve  which  was  measured  in  an 
experiment.  Let  us  assume  that  we  have  a  uniaxial  stress-strain  curve,  as  shown  in 
Figure  1.  Stress  is  written  by  the  following  equation: 

cr  =  (\  -  d  )Es  (22) 

Therefore,  the  damage  parameter  d  is  expressed  by  the  stress  and  strain,  as: 

d  =  1-—  (23) 

Es 

The  value  of  the  effective  stress-like  parameter  f  can  be  written  in  terms  of  the 
uniaxial  stress  and  strain,  as: 


T  = 


y[a°£  = 


(24) 


Therefore,  once  the  uniaxial  stress-strain  curve  is  given,  relationship  between  the 
effective  stress-like  and  the  damage  parameter  can  be  obtained.  In  present  investigation, 
a  uniaxial  stress-strain  curve  is  approximated  by  a  series  of  piecewise  straight  lines,  as 
depicted  in  Figure  2.  In  Figure  2,  points  on  the  stress-strain  curve  (a j,£X),  (<j2 ,s2), 
(cr3,£3),  •••,  (a; (<j/+1  , si+l ) ,  •••  connect  straight  line  segments.  From  a  set  of 
data  (crj,^),  (<r2,e2),  (<r3,e3) ,  •••  ,  (cr,-,ef) ,  (eri+1,si+1) ,  •••  ,  we  can  compute 
(r, , ,  (f2,d2),  (r3 ,d3),  •••,  (f, , dj ) ,  •••  by  using  equations  (23)  and 

(24).  Thus,  the  relationship  between  f  and  d  is  also  approximated  by  a  series  of 
piecewise  straight  lines.  The  function  H(f,d)  which  governs  the  damage  evolution  law 
of  equation  (6)  is  approximated  by: 

h(t, d)  =  4- « -^4 — zr~  (dj  <  d  <  dj+i,Tj  <  t  <  ri+i)  (25) 

r  Ti+ 1  -  Ti 

The  function  H(f,d)  whose  value  is  determined  by  equation  (25)  is  used  in  the 
stress-strain  relationship  of  equation  (10). 
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Slope:  £  (Young's  modulus) 


Fig.  1  Uniaxial  stress-strain  curve 


Slope:  E  (Young's  modulus) 


Fig.  2  Piecewise  linear  approximation  for  uniaxial  stress-strain  curve 


1.4  Determining  damage  evolution  law  for  the  separate  Isotropic/deviatoric 
damage  model 

In  this  section,  the  damage  evolution  law  for  the  separate  isotropic/deviatoric  damage 
model  is  dealt  with.  As  seen  in  equations  (17)  and  (18),  there  are  two  different  damage 
evolution  laws  and  two  damage  variables  Hv  (fv ,  dv )  and  HD (fD ,  dD ) .  In  order  to  fully 
determine  Hv  (fv ,  dv )  and  HD  (rD , dD ) ,  we  need  at  least  two  sets  of  stress-strain  curves. 
When  only  one  stress-strain  curve  is  available,  we  need  to  introduce  at  least  an 
additional  condition.  Procedures  to  determine  the  damage  evolution  laws  from  only  one 
set  of  uniaxial  stress-strain  curve  are  described  in  this  section. 

We  assume  that  a  uniaxial  stress-strain  curve,  as  shown  in  Figure  1,  is  available.  As  an 
additional  condition  to  uniquely  determine  the  damage  evolution  law,  we  postulate  that 
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the  ratio  ( dv  jdD  )  between  the  dilatational  and  deviatoric  damage  parameters  remains  to 
be  the  same  during  the  uniaxial  deformation.  We  write: 

d  q  =  ccdy  (26) 

Here,  a  is  a  positive  constant  ( 0  <  a  <  oo ).  When  a  equals  zero,  only  the  dilatational 
damage  is  present.  When  a  =  1 ,  the  magnitudes  of  dilatational  and  deviatoric  parts 
equal  each  other  and,  therefore,  this  case  is  very  similar  to  that  of  the  isotropic  damage. 
When  a  is  infinitely  large,  the  dilatational  damage  parameter  dv  is  zero.  Therefore, 
material  undergoes  deviatoric  damage  only. 

We  assume  a  uniaxial  deformation  in  x,  direction.  We  can  write  the  following 
statements. 


cjji  *  0  (unknown),  o2 2  =  <J33  =  <r12  =  a23  =  cr31  =  0  (known) 

£|  |  i=-  0,  £'i 2  =  £'23  =  £31  =  0  (known)  s22  —  s33  (unknown) 


(27) 


Therefore,  by  substituting  equations  (26)  and  (27)  in  equation  (19),  we  have: 


°il  -  +  £22  +  ^33)+  ~(l_ adv)/j(2sn  -£22  “^33) 

2 

°22  =  (l  -  dy  )K{£\  1  +  £'22  +  f33  )  +  ^  (l  _  adV  )t-l(~£22  “f33  _£il) 

2 

°il  =  (1_^vK(fll  +£22  +^33)+—  (\-adv)ll{2.£33  -£\1  -£ 22 ) 


(28) 


From  the  second  and  the  third  of  equation  (28)  and  e2 2  =  s33 ,  one  can  obtain: 


0  -  (l  -  dv  )K(&'}  1  +£22  +  £33 )+  “0  “  ady  )m(£]  1  -£22) 


(29) 


The  first  of  equation  (28)  and  equation  (29)  lead  to: 


£22  ~ 


-Qii 

2(l  -  ady  )// 


+  £li 


(30) 
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By  substituting  equation  (30)  in  the  first  of  equation  (28),  we  can  establish  a 
relationship  between  crn  and  , ,  as: 

?>{^-ady^-dy)jjK£ii  -  ji(l-ady)//  +  (l-dy)/rj(Tii  =0  (31) 


By  rearranging  equation  (31),  we  can  establish: 
a/j  +  3K  cjjj 


ad2+)^^L^lL.(l  +  a)U 
1  9  juK  su  v  /( 


f  gjj  ^ 
Esu 


1- 


=  0 


(32) 


Therefore,  for  a  given  set  of  variables  an ,  £■, ,  and  «  ,  we  can  solve  for  dv . 
Two  special  cases  ( a  =  0  and  a  =  1 )  are  presented  first.  When  a  =  0 ,  we  have: 


O  dy  + 


0  ■  //  +  3K  a | 
9/jK 


f-(i  +  o)Uv 


fdy  + 

fi-  <T>>  1 

J 

v  1  y 

=  0 


and. 


e?v  —  1  ■ 


3a: 


q~ii 
3^ii  - 


o-il 


(a  =  0) 


(33) 


(34) 


When  a  =  1 ,  we  have: 

^  +  crn 


4  + 


9//AT  £-n 


Wy  + 

fi- a" ) 

1 

{  Een) 

=  0 


and. 


dv  —  1 — 


°il 
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(a  =  0) 


(35) 


(36) 


In  is  noted  that  there  are  two  solutions  that  satisfy  equation  (35)  and  they  are 

dv  =  1  -  -^LL  and  dv  =  1 .  The  second  one  takes  a  constant  value  and,  therefore,  is  not  a 
Es  ii 

feasible  solution.  For  more  general  case  (a*  0,1 ),  we  solve  equation  (32)  for  dv ,  by 
letting: 

ccdy  +  Ady  +  3  —  0  (37) 
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A  =  m±lL?n-(Ua)  and  B  =  i 
9  juK  eu 


(38) 


^11 

Esxx 


Therefore,  the  value  dv  is  determined  to  be: 


d\r  — 


-A+Va2  - 


4  aB 


2a 


(39) 


In  equation  (39),  the  sign  associated  with  Va2  -4aB  needs  to  be  determined.  To 

determine  the  sign,  we  check  a  special  case  ( a  =  1 ).  By  letting  a  =  1 ,  in  equation  (39), 
we  have: 


dy 


1  1 

[ 

(  <T\\ 

\ 

-2 

|+  1 

2  1 

[ 

,Es\\ 

J 

1  Esn 

(40) 


When  we  take  the  positive  sign,  equation  (40)  results  in: 


1  I 

f 

A 

-2 

U  ""  ) 

2  j 

l^n 

J 

1  Esn 

(41) 


This  is  a  constant  value  ( dv  =  1 )  and,  therefore,  is  inappropriate.  When  the  negative  sign 
is  assumed,  we  obtain: 


dy 


1  I 

f 

A 

-2 

an 

2  j 

r 

i^n 

J 

Esn 

(42) 


This  result  is  the  same  as  the  case  of  isotropic  damage  material  (equation  (23)). 
Therefore,  we  choose  the  negative  sign  in  equation  (39). 


dy  — 


(43) 


Thus,  the  deviatoric  damage  parameter  dD  is  also  determined  to  be: 
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(44) 


Once  uniaxial  stress-strain  curve  as  depicted  in  Figure  1  is  given,  we  collect  a  series  of 
points  (ax,sx),  (<7 2 , £'2 ) ,  (<73,£-3),  •••,  (cr, , £, ) ,  (aM,£M) ,  •  ••,  on  the  curve,  as 

shown  in  Figure  2.  Damage  parameters  ( dv\.  and  dD |.)  corresponding  to  the  stress 

and  strain  (a are  obtained.  Thus,  by  using  equations  (27)  and  (28),  unknown 
strains  £22  and  e33  ( £-22  =  s33 )  are  derived.  Therefore,  we  can  compute  the  effective 

stress-like  terms  tv  (=  4Kskk )  and  rD(=  .  Thus,  we  have  a  series  of  data 

(dV\vdD\vfy\vfD\x)  ,  (dy\2,dD\2,fV\2,fD\2  )  .  H  ,dD\3,Tv\3,TD\3)  ,  •••  , 

(dv\.,dD\rTV\.,TD\.),  (dy |.+1  ,dD | f+1 , fy | ,+1  ,rD|.+|),  ....  The  scalar  functions,  Hv (rv , dv ) 
and  HD(fD,dD )  are  derived  to  be : 
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2.  Numerical  Implementation  of  the  Damage  Constitutive  Laws  in  an 
In-House  Finite  Element  Program 

The  isotropic  and  the  separate  isotropic/deviatoric  damage  models  were  implemented  in 
an  in-house  finite  element  program.  In  this  chapter,  the  numerical  implementations  of 
the  damage  constitutive  laws  are  described.  The  in-house  finite  element  program  that  is 
used  in  present  investigation  is  based  on  the  s-version  finite  element  method  (s-FEM). 
In  this  chapter,  (i)  the  in-house  finite  element  program  based  on  s-FEM,  (ii)  the 
implementation  of  nonlinear  constitutive  laws  in  the  s-FEM  program,  (iii)  the 
implementation  of  isotropic  damage  constitutive  model  and  (iv)  the  implementation  of 
isotropic/deviatoric  damage  model  are  discussed. 

2.1  The  s-version  finite  element  method  (s-FEM) 

In  present  investigation,  the  s-FEM  program  that  has  been  developed  in  Kagoshima 
university  was  adopted  as  the  base  program.  As  described  in  Okada  et  al.  [3]  and  Tanaka 
et  al.  [4],  the  s-FEM  is  the  most  useful  in  the  analyses  of  composite  material.  However, 
the  code  can  be  used  as  an  ordinary  finite  element  program.  In  this  section,  the  s-FEM 
for  nonlinear  analysis  is  briefly  described. 

The  s-FEM  was  first  proposed  by  Fish  [5]  and  has  been  applied  to  various  engineering 
problems  [5, 6, 7, 8].  In  the  s-FEM,  two  types  of  finite  element  models  are  used.  One  is 
called  “global  finite  element  model”  or  “global  model”  that  covers  the  whole  analysis 
region.  The  other  is  called  “local  finite  element  model”  or  “local  model”  that  is 
superposed  on  the  global  model,  as  depicted  in  Figure  3.  Typically,  the  local  model  has 
a  finer  spatial  resolution  than  the  global  model  and  is  typically  placed  at  the  portion  of 
stress  concentration,  such  as  crack  (see  Okada,  Endoh  and  Kikuchi  [7]),  and  second 
phase  material  and  its  vicinity  (see  Okada  et  al.  [8]  and  Okada  et  al.  [9]).  In  Okada  et  al. 
[3],  an  s-FEM  formulation  that  allowed  the  local  models  to  multiply  overlap,  as  shown 
in  Figure  4. 

Let  us  denote  the  regions  of  global  model  and  local  models  to  be  nG  and  Q.Ll 
(/  =  1,2, 3, •••,«),  respectively.  Here,  as  done  in  Okada  et  al.  [3],  there  are  many  local 

models  that  are  superposed  on  the  global  model  and  they  are  allowed  to  multiply 
overlap  each  other.  Displacements  are  expressed  by  the  shape  functions  of  finite 
elements  of  global  and  local  models  independently.  This  means  that  we  do  not  need  to 
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maintain  any  relationships  between  the  locations  of  nodes  and  elements  of  the  global 
and  local  models.  Therefore,  it  is  quite  tractable  to  enhance  the  spatial  resolution  of 
analysis  model  locally,  by  superposing  the  local  models  on  the  global  model.  It  is  noted 
that  if  there  were  no  overlaid  local  models,  the  program  performs  an  ordinary  finite 
element  analysis. 

When  the  damage  constitutive  law  is  adopted  in  an  analysis,  we  perform  an  incremental 
analysis,  just  like  the  case  of  elastoplasticity  (see  Okada  et  al.  [9]).  We  briefly  describe 
the  formulations  of  present  s-FEM. 

We  denote  the  displacement  increments  to  be  A uf  and  Auf  (/  =  1,2,3, •••,«)  that  are 
based  on  the  shape  functions  of  finite  elements  of  the  global  and  the  local  models.  The 
displacement  functions  are  superposed  when  the  global  and  the  local  mesh  regions 
overlap  each  other.  For  example,  for  a  region,  where  the  global  and  the  local  mesh 
regions  ( Q.Lp ,  Q.Lq  and  C±Ll )  overlap,  we  write  the  displacement  increments  Am,  ,  as: 

Am,=  Auf  +  Am  9’  +  Aufq  +  Auf1  (50) 


At  a  point,  where  any  local  models  do  not  overlap  with  the  global  model,  the 
displacements  w,  equal  uf  ( ut  =  uf ).  The  continuities  of  the  displacements  are 
enforced  by  letting  the  functions  uf  ,  uf2,  uf 3,  ••• ,  ufn  be  zero  at  the  outer 
boundaries  rLi  of  the  local  model  regions  Q.Ll .  The  statement  of  principle  of  virtual 
work  in  an  incremental  formulation  is  written  to  be: 


dASu 
dx , 


i  Damage  SAu^ 

W  ~dfff 


dQ6’  =  fa  SAU'At^rf 


(51) 


where  the  variations  of  the  displacements  are  defined  in  the  same  manner  as  the 
displacements,  r f  denotes  the  boundary  where  the  tractions  r“  are  prescribed. 

D*?*  are  t^e  forth  order  tensor  representing  material’s  stress-strain  relationship 

(incremental  constitutive  law).  We  then  substitute  the  displacements  and  their  variations 
in  the  statement  of  principle  of  virtual  work.  After  some  algebraic  manipulations,  we 
have: 
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dASu 


dAu 


JQ 


jj  Damage 

dx,  w*  dx 


k  dOG 


U< 


q=l 


0A  Su? 
dx ; 


p.  Damage  ^Au^ 

W  —fac 


Lq 


-dQ.1*  =  |r6-  Su^At^vf1 


In4’ ' 


dASu , 


Lp 


dx ; 


i  Damage 


ijki 


dAuLk 

dxo 


-dO.Lp  +|, 


hLp 


Damage  5Ai^  dQ Lp 
8X  :  <JU  8X , 


n  dASu^ 

■  E  \QLP-Lq— 

q= 1  j 


z  ^  Damage 


8Au^q  ,  , 

— — ciql,,_Zy/  =o 

dxt 


(p  =  1,2,3,  •••,«) 


(52) 


(53) 


It  is  noted  here  that  the  right  hand  side  of  equation  (53)  equals  zero,  since  the  local 
model  regions  do  not  intersect  with  the  outer  boundary  of  the  structure.  nLp~Lq 
(  p,q  =  1,2,3,- ••,«;  p *  q  )  denote  regions  sheared  by  fiLp  and  Q Lq  (i.e., 
nLp~Lq  =qLP  n q/-</  j  when  o  =  nLp  r,nLq ,  the  related  terms  in  equation  (53)  disappear. 
After  the  global  and  local  regions  are  discretized  by  appropriate  finite  elements,  we  can 
write  a  matrix  formulation,  as: 


A fG 

0 

0 

>  = 

G  G—Ll  G-L2 

j^Ll-G  kL\  g  L1-L2 

kL2-G  L2-L1  kL2 

G-Ln 

g  Ll-Ln 

L2—Ln 

A  uG 

A uLl 
Au 12 

0 

Ln-G  Ln-Ll  Ln—L2 

..  KLn 

Au Ln 

(54) 


where  A uG  and  A uLp  are  the  column  vectors  of  unknown  nodal  displacements  of  the 
global  and  the  local  meshes.  Matrices  kg  ,  KLp ,  kg~Lp  ,  KLp~Lq  arise  from  the 

integrals  in  the  left  hand  sides  of  equations  (52)  and  (53).  Since,  kg~Lp  =  (kLp~gJ  and 

KLP-Lq  _  (k  J  ;  the  global  stiffness  matrix  in  the  right  hand  side  of  equation  (54)  is 

symmetric.  However,  the  global  stiffness  matrix  is  not  banded,  unlike  an  ordinary  finite 
element  method.  A fG  is  the  consistent  nodal  force  vector  arising  from  the  right  hand 
side  of  equation  (52). 

The  coefficient  matrix  of  equation  (54)  is  symmetric.  However,  the  coefficient  matrix 
looses  its  band  structure,  as  there  are  many  off-diagonal  components.  In  such  a  case,  use 
of  iterative  equation  solver  is  advantageous  over  direct  solvers.  Thus,  we  adopted  an 
Incomplete  Cholesky-decomposition-preconditioned  conjugate  gradient  (ICCG)  method 
[10]. 
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Figure  3  A  schematic  illustration  for  s-FEM  (s-version  finite  element 
method)  modeling.  A  crack  problem  is  presented  as  an  example. 


(a)  s-FEM  discretization  (b)  Global  and  local  model  regions 

Figure  4  The  s-FEM  model  with  multiple  local  model  regions  that  overlap  each  other. 


2.2  Incremental  algorithm 

The  constitutive  equations  were  implemented  in  the  s-FEM  computer  program.  In 
present  analysis,  we  did  not  implement  equilibrium  iteration  scheme.  The  reason  why 
we  did  not  do  so  was  that  computations  for  nodal  reaction  forces  became  sometimes 
erroneous.  Hence,  equilibrium  iteration  algorithms  such  as  Newton-Raphson  scheme 
may  not  converge  or  may  converge  to  wrong  solutions.  The  displacements  w,(x), 

stresses  <r(/  (.v)  and  strains  £ij(x)  at  a  point  are  expressed  by  the  sum  of  their 
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increments  A m/(x),  A cr/(x)  and  A*r/(x),as: 


Uj  (x)  =  X  Am/  (x)  ,  & ij  (x)  =  x  Act/  (x)  and  (x)  =  X  A £■:  (x)  (55) 

/  I  1 

All  the  other  deformation  parameters  are  treated  in  the  same  fashion. 

2.3  Evaluation  of  relationship  between  the  effective  stress-like  parameters  and  the 
damage  parameters 

The  procedures  that  are  described  in  sections  1.3  and  1.4  are  implemented  in  the  s-FEM 
computer  program.  Therefore,  an  analyst  only  needs  to  specify  a  one  dimensional 
stress-strain  relationship,  Young’s  modulus  and  Poisson’s  ratio.  The  stress-strain 
relationship  is  given  by  specifying  a  series  of  data  points  on  the  curve  [  (ox ,  sx ) ,  (o2  ,s2), 
(o3,s3),  •••,  (o/ , 6'j ) ,  (oM , ) ,  •••].  The  first  data  point  (cr, , £\ )  must  be  on  the 
line  of  elastic  slope,  i.e.;  o1=Ee1,  where  E  is  the  Young’s  modulus. 

Thus,  the  analyst  does  not  have  to  deal  with  the  tedious  procedures  of  sections  1.3  and 

1 .4  and  the  input  data  structure  is  very  analogous  to  that  of  isotropic  elastoplasticity. 

2.4  Some  other  issues  associated  with  the  damage  constitutive  law-initiation  of 
nonlinear  deformation 

Material  initially  undergoes  elastic  deformation  and  then  follows  nonlinear  deformation 
path.  That  is  analogous  to  the  case  of  elastoplasticity.  In  this  section,  algorithms  when 
the  nonlinear  deformation  initiates  at  a  point  during  an  increment  are  discussed.  First, 
we  deal  with  the  case  of  isotropic  damage. 

Let  us  assume  that  uniaxial  stress-strain  curve  is  given  as  shown  in  Figure  5.  The 
stress-strain  relationship  upto  (o} ,  sx )  is  linear  (no  damage)  and  the  nonlinear 
deformation  initiates  at  (ox,£X).  Then,  by  following  the  procedures  in  section  1.3,  the 
sets  of  effective  stress  and  damage  parameter  corresponding  to  ( <j2,£2 ),  {o3,e3),  •••, 
(cr, , £, ) ,  (cr,+| , £[+\ ) ,  •••  are  determined  to  be  (r,, <:/,),  {f2,d2),  (f3,d3),  •••,  (r(- , dt ) , 

(fM,dM),  •••  .  We  assume  that  the  state  of  stresses  cr/  and  strains  efj  at  the 
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beginning  of  /  th  incremental  step  satisfies: 


^ijk(.£ij£ki  <  rl 


I  T 


(56) 


It  is  also  assumed  that  during  the  I  th  incremental  step,  material  damage  starts  taking 
place.  Therefore,  we  can  write: 


Cijkt 


(  Sij  ^£ij  \£k(  "*■  ^£k(  ) 


>n 


(57) 


Thus,  the  I  th  incremental  step  can  be  split  into  elastic  and  nonlinear  part.  We  let  the 
their  fractions  to  be  a  and  1  -  a  .  Hence,  the  following  equation  holds: 


Cijkt ' 


|  £ jj  +  CxASjj  k,  +  CcASfcp  T\ 


(58) 


From  equation  (58),  we  can  determine  the  value  of  a  . 

0  =  Cijkl  ( £ij  +  a^£ij  \£ki  +  a^-£ki)~^\ 

=  ifijM£ij£U  -  )+  2aCijke£ijA£k(  +  a2cijkeA£ijAsk£  (59) 

=  C  +  2,ccB  +  cc  A 

where  C<0  and  A > 0 ,  since  the  damage  state  has  not  yet  been  reached  and  the 
quadratic  form  of  CijM  is  positive  definite,  a  can  be  determined  to  be: 


a  ■ 


■  2 B  +  V4B^4AC  -  B  ±  Vfi2  -  AC 


2A 


(60) 


Since  C<0  and  A>0,  AC<0.  We  take  the  positive  sign,  in  order  for  a  to  be  a 
positive  number.  Therefore,  we  have: 


a  - 


-B+Vfi2 -AC 


(61) 
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Thus,  the  increments  Aa \  of  stresses  during  the  I  th  incremental  step  can  be  written 


to  be: 


Acri j  -  Cijkl 


(aAsjj )+  (l - aj]Cijkl Ae 


+  ocC^  ■■  A  f 


Jcr  [(  +  aCkipqA£pq  ]asI(  1  (62) 


For  the  separate  Isotropic/deviatoric  damage  model,  the  similar  procedures  can  be 
established.  They  are  described  below. 

First,  we  assume  that  a  uniaxial  stress-strain  curve  as  depicted  in  Figure  5  is  given.  We 
extract  a  series  of  data  on  effective  stresses  and  damage  parameters  {dv\vdD\vfv\vzD\x), 

[dv\2,dD\2,fy\2,fD\2)  ,  (dV\3,dD\3,¥v\3,TD\3)  ,  •••  ,  (dv\.,dD\.,TV\rTD\.)  , 

(dv\i+vdD\.+vTV\.+vTD\.+l) ,  •••  from  the  stress  and  strain  data  (al,el) ,  (<r2 . £'2 ) , 
(0-3 , £3 ) ,  •••,  (a (<t,+| , el+\ ) ,  •••.  We  assume  that  the  states  of  stresses  and 
strains  ( ajj ,  slq )  and  ( a-j ,  s-j )  at  the  beginning  of  I  th  and  J  th  incremental  step 
satisfy: 


^^£ij  £ij  <  \TV  j  ) 

(63) 

K{dkf  <fe1)2 

(64) 

Then,  it  is  also  assumed  that,  during  the  I  th  and  J  th  incremental  step,  deviatoric  and 
dilatational  damage  start  taking  place.  Thus,  we  can  write: 


+A £ij\sij  +  Ae’d ) >  (z>  | L )  (65) 

K{skk  +  A6"kk  J  > 


21 


(66) 


The  I  th  and  J  th  increments  are  split  into  two  portions  that  are  elastic  and  nonlinear 
parts.  The  fractions  of  the  elastic  parts  are  specified  to  be  a7  and  a7  .  Thus,  we  can 
write: 


2/4/  +alA£,ij  jelj  +  a1  As\J )=  (tD\x  f 


KVkk  +  a  A £kk 


Then,  a1  and  aJ  can  be  determined.  From  equation  (67),  we  obtain: 


(a1  J2 Ae[j  Ae\j  +2a7 s\j As'j  +  s'J  £'j  -^-  =  0 


(67) 

(68) 


(69) 


By  letting  A°  =  A e-j A s\j  ,  B°  =  £\J A£\j  and  C  =  e'J £\j 


we  arrive  at: 


a1  = 


-2 Bd  ± 


UnDJ  - 


4ADCD  -Bd± 


Ml 


adcd 


2  A 


D 


iD 


(70) 


In  equation  (70),  it  is  appropriate  to  choose  the  positive  sign.  Therefore,  we  have: 


a 1  = 


-  bd  +  ^(bd  J  -adcd 


i  D 


(71) 


From  equation  (68),  we  obtain: 


(aJf( A4f +2«j4a4  + 


kkf- 


(Mi)2 


K 


=  0 


(72) 


Thus,  equation  (72)  is  solved  for  aJ  by  letting  Av  =[ A£Jkkf  ,  Bv  =  4A4  ar>d 

(jf  (fv|,)2 

\£kk )  v 


cv  = 


and  we  have: 
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(73) 


In  equation  (73),  the  appropriate  sign  associated  with  -AVCV  has  already  been 

chosen.  The  increments  of  the  deviatoric  and  hydrostatic  stresses  during  the  /  th  and 
the  J  th  incremental  steps  are  written  to  be: 


A<4  =2(1  -d^’j 

rD 

=(i-^)*:a4 -d-»)^fejv)(g4fA4 


Slope:  E  (Young's  modulus) 


£ 


Figure  5  An  illustration  of  stress-strain  curve  and  its  transition  from  the  linear 
(elastic)  to  nonlinear  (damage)  deformation. 


23 


3.  Numerical  demonstrations  of  s-FEM  computer  program  with  the 
damage  constitutive  laws 


In  this  chapter,  the  s-FEM  program  that  is  developed  during  this  course  of  study  is 
demonstrated.  Example  problems  that  contain  one  or  two  particles  in  a  unit  cell  region 
are  solved. 

3.1  Stress-strain  relationship 

Stress-strain  curve  was  taken  from  Kwon  and  Liu  [11]  and  is  shown  in  Figure  6.  We 
then  extract  a  series  of  data  points  (gx,£\),  (<t2 , £2 ) ,  (<x3 ,  <?3 ) ,  •  •  •  ,  (a 7,e7),  as  also 
shown  in  Figure  6.  As  an  example,  the  relationship  between  the  uniaxial  strain  and  the 
damage  parameter  for  the  case  of  isotropic  damage  is  depicted  in  Figure  7.  Thus,  the 
uniaxial  stress-strain  curve  is  reconstructed  by  a  series  of  linear  approximations,  as 
shown  in  Figure  8. 
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Fig.  6.  Stress-strain  curve  of  a  uniform  composite  specimen. 


Figure  6  Uniaxial  stress-strain  curve  of  polymeric  particulate  composite  material 
that  was  given  in  Kwon  and  Liu  [11] 


24 


0  0.05  0.1  0.15  0.2  0.25 

Strain 


Figure  7  Relationship  between  the  strain  and  the  damage  parameter  for  the  isotropic 
damage  constitutive  model,  which  was  extracted  from  the  stress-strain  curve  of 
Figure  6. 


Figure  8  Uniaxial  stress-strain  curve  that  was  reconstructed  after  the  relationship 
between  the  effective  stress-like  and  the  damage  parameters  is  established  by  the 
proposed  procedures. 


3.2  Example  problem-isotropic  damage  model 
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In  this  section  a  simple  example  problem  that  demonstrates  the  isotropic  damage 
constitutive  model  is  presented.  A  block  containing  two  hard  particles  subject  to  tension, 
as  depicted  in  Figure  9,  is  solved.  To  model  it,  we  used  the  s-FEM  approach.  The 
particle  and  its  vicinity  were  discretized  by  local  model  and  the  block  as  whole  was 
modeled  by  the  global  one.  They  are  shown  in  Figure  10. 


In  Figure  11-14,  the  distributions  of  stress  and  the  evolution  of  damage  parameter  at  a 
section  cutting  through  the  center  of  the  particle  are  shown.  It  is  seen  that  the  material 
damage  develops  at  the  top  and  bottom  of  the  particles.  The  damage  region  connects  the 
particles. 


Tension  a 


Figure  9  A  cubical  block  containing  two  particles,  subject  to  tensioin. 


(a)  Global  finite  element  model  (b)  Local  finite  element  model 

Figure  10  The  global  and  local  finite  element  models  for  the  s-FEM  analysis. 
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Figure  1 1  Tensile  stress  distribution 
at  the  tensile  strain  6.32x1  O'2 


Figure  12  Distribution  of  the  damage 
parameter  at  the  tensile  strain  6.32x1  O'2 
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Figure  13  Tensile  stress  distribution  Figure  14  Distribution  of  the  damage 
at  the  tensile  strain  lO.OxlO"2  parameter  at  the  tensile  strain  10.0x1  O'2 


3.3  Example  problem-separate  Isotropic/deviatoric  damage  model 

The  separate  isotropic/deviatoric  damage  model  is  also  tested  for  the  same  problem.  The 
constants  a  which  determines  the  ratio  of  contributions  of  hydrostatic  pressure  and 
deviatoric  stress  was  varied  (0.1,  1.0,  10.0).  When  a  =  0.1,  the  contribution  of  the 
hydrostatic  pressure  dominates  that  of  the  deviatoric  stresses.  When  a  =  1.0,  the 
material  response  should  be  very  similar  to  that  of  the  isotropic  damage  case.  When 
a  =  10.0 ,  the  deviatoric  stress  components  should  have  the  major  contribution. 

In  Figures  15,  16  and  17,  the  distributions  of  stress  and  the  damage  parameters  in  a 
section  that  cuts  through  the  center  of  the  particles  are  depicted  for  a  =  0.1,  a  =  0.99 
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and  a  =  10.0,  respectively.  The  reason  why  we  did  not  let  a  be  1.0  was  that  a  =  1.0 
requires  a  special  program  implementation  in  the  procedures  which  are  described  in 
section  1.4. 

In  the  case  of  a  =  0.1 ,  it  is  seen  that  the  magnitude  of  dilatational  damage  parameter  is 
much  larger  than  that  of  the  deviatoric  one.  Therefore,  as  we  intended,  the  dilatational 
damage  is  the  major  damage  mode.  When  a  =  0.99 ,  the  magnitudes  of  the  dilatational 
and  deviatoric  damage  parameters  are  similar.  For  a  =  10.0 ,  the  value  of  the  deviatoric 
damage  parameter  is  much  larger  than  that  of  the  dilatational  one.  In  this  case  also,  the 
material  behaves  as  we  intended. 
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Figure  15  The  distributions  of  tensile  stress  and  the  damage  parameters  in  the  section 
cutting  through  two  particles  when  the  separate  dilatational/deviatoric  damage 
constitutive  model  with  a  =  0.1  is  used. 
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Figure  16  The  distributions  of  tensile  stress  and  the  damage  parameters  in  the  section 
cutting  through  two  particles  when  the  separate  dilatational/deviatoric  damage 
constitutive  model  with  a  =  0.99  is  used. 
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Figure  17  The  distributions  of  tensile  stress  and  the  damage  parameters  in  the  section 
cutting  through  two  particles  when  the  separate  dilatational/deviatoric  damage 
constitutive  model  with  a  =  10.0  is  used. 


3.4  Summary-Numerical  demonstrations  on  the  damage  constitutive  models  and 
the  s-FEM  computer  program 

The  damage  constitutive  models  are  tested  in  this  chapter.  We  found  that  material 
damages  accumulate  at  the  locations  of  the  stress  concentration.  The  damage  regions 
connect  between  the  particles. 

The  separate  isotropic/deviatoric  damage  model  could  control  the  damage  mode.  The 
intended  damage  mode  (dilatational  or  deviatoric)  dominates  the  other. 
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4.  Investigating  the  Effects  of  Micro-Structure  on  the  Tensile  Behavior 
of  a  Solid  Propellant 

Having  developed  the  damage  models  and  the  s-FEM  computer  program,  we  carried  out 
a  series  of  analyses  on  the  tensile  behavior  of  a  block  made  of  polymeric  particulate 
composite  material.  Thought  the  particles  distribute  evenly  in  a  macroscopic  sense,  they 
have  some  distributions  at  a  microscopic  level.  In  order  to  investigate  the  effects  of  the 
locally  uneven  distributions  of  the  particles,  we  assume  that  the  variations  in  material 
stiffness  at  a  local  level.  Experimental  evidences  showed  that  when  representative 
volume  of  larger  than  2.0  mm  by  2.0  mm  the  material  could  be  treated  as  a 
homogeneous  material.  This  means  that  if  one  took  the  representative  volume  equal  to 
or  larger  than  2.0  mm  by  2.0  mm  the  volume  fraction  of  particles  is  approximately  the 
same  throughout  the  block. 

The  reference  stress-strain  curve  was  taken  from  Kwon  and  Liu  [11].  Magnification 
factors  are  multiplied  and  six  different  stress-strain  curves  (level  -2.5,  -1.5,  -0.5,  0.5,  1.5, 
2.5)  were  generated  as  shown  in  Figure  18.  Then,  by  using  a  random  number  generator 
with  Gaussian  distribution,  a  series  of  random  numbers  are  generated  and  are  assigned 
to  finite  elements.  Material  data  is  assigned  to  each  finite  element  according  to  the  value 
of  assigned  random  number.  For  example,  when  the  value  of  assigned  random  number 
is  between  -p  and  -p/2 ,  the  stress-strain  curve  of  level  -1.5  is  assigned  to  the 
element.  The  relationships  between  the  assigned  random  number  and  the  level  of 
stress-strain  curve  to  be  used  are  shown  in  Figure  19  and  Table  1. 


Table  1  The  relationship  between  the  assigned  random  number  and  the  level  of 

stress-strain  curve  to  be  adopted. 


Random 
number  V 

V  <-p 

-p  <V  <  -0.5  p 

-0.5/9  <  V  <  0 

o<y  <0.5  p 

0.5 p  <  V  <  p 

p<V 

Level 

-2.5 

-1.5 

-0.5 

0.5 

1.5 

2.5 

Magnification 

Factor 

0.75 

0.85 

0.95 

1.05 

1.15 

1.25 

Two  different  analysis  models  were  generated  by  specifying  different  initial  values  to 
the  random  number  generator.  They  are  different  but  are  statistically  the  same.  The 
models  are  shown  in  Figure  20.  The  finite  element  model  is  shown  in  Figure  21.  There 
are  20000  finite  elements.  The  size  of  finite  element  is  less  than  0.2  mm.  Therefore, 
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there  are  at  least  100  elements  inside  the  2.0  mm  by  2.0  mm  representative  area. 
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Reference  curve 


Fig.  6.  Stress-strain  curve  of  a  uniform  composite  specimen. 
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Figure  1 8  Six  level  stress-strain  relationships  that  were  postulated  in  present 
investigation. 


Figure  19  Relationship  between  the  levels  of  the  postulated  stress-strain  curves 
and  the  value  of  generated  random  number. 
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Figure  20  Two  different  analysis  models  that  have  the  statically  the  same 
distributions  of  material  stiffness. 


(a)  Tension  block 


(b)  With  a  crack 


Figure  21  The  finite  element  analysis  models 


4.1  The  behavior  of  a  tensile  block-isotropic  damage  constitutive  model 


In  this  section,  the  tensile  behavior  of  a  block  that  was  made  of  the  isotropic  continuum 
damage  solid  is  presented.  Stress-strain  curve  at  a  material  point  was  set  by  the 
procedures  that  are  presented  in  the  previous  section.  The  block  was  subject  to  15%  of 
tensile  strain. 
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In  Figure  22,  the  distributions  of  stress  in  the  tensile  direction  when  the  block  is  subject 
to  0.3%  of  the  strain  are  depicted.  In  this  step,  all  the  material  points  in  the  blocks  have 
not  undergone  any  material  damage.  Stress  distributions  show  that  the  regions  of  low 
and  high  stresses  run  in  the  direction  of  loading. 

In  Figure  23,  the  distributions  of  the  stress  when  the  block  is  subject  to  10%  of  tensile 
strain  are  shown.  The  trends  in  the  figures  are  the  same  as  those  in  Figure  22.  In  Figure 
24,  the  distributions  of  the  isotropic  damage  parameter  are  presented.  Figure  24  shows 
that  major  trends  in  two  different  statistically  the  same  blocks  are  the  same.  The  regions 
of  relatively  sever  damage  run  in  about  45  degree  from  the  axis  of  loading.  However,  it 
is  seen  that  the  regions  of  very  sever  damage  run  perpendicular  to  the  tensile  axis. 

In  Figure  25,  the  distributions  of  the  stress  when  the  block  is  subject  to  15%  of  strain 
are  depicted.  The  major  trends  are  the  same  as  Figures  21  and  23.  The  regions  of  high 
and  low  stresses  run  vertically.  The  distributions  of  the  damage  parameter  are  shown  in 
Figure  26.  The  major  trends  are  the  same  between  two  cases  that  are  statistically  the 
same.  It  is  seen  that  many  regions  of  sever  damage  run  in  short  distances,  typically 
equal  or  shorter  than  1  mm,  perpendicular  to  the  direction  of  loading.  Then,  they  link 
each  other  making  larger  damaged  regions.  The  larger  damaged  regions  run  in  various 
directions. 
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Figure  22  Tensile  stress  distributions  when  the  blocks  are  subject  to  0.3%  of  strain 

(elastic  stage). 
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Figure  23  Tensile  stress  distributions  when  the  blocks  are  subject  to  10%  of  strain. 


(a)  Model  #1 


(b)  Model  #2 


Figure  24  The  distributions  of  the  damage  parameter  when  the  blocks  are  subject  to 

10%  of  strain. 
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Figure  25  Tensile  stress  distributions  when  the  blocks  are  subject  to  15%  of  strain. 
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Figure  26  The  distributions  of  the  damage  parameter  when  the  blocks  are  subject  to 

15%  of  strain. 


4.2  The  behavior  of  a  tensile  block-  separate  Isotropic/deviatoric  damage  model 


In  this  section,  the  tensile  behavior  of  the  block  that  was  analyzed  by  the  separate 
isotropic/deviatoric  damage  model  is  presented.  The  parameter  a  which  characterizes 
the  ratio  of  contribution  of  the  volumetric  and  deviatoric  damage  modes  was  set  to  be 
0.1,  0.99  or  10.0.  These  values  correspond  to  the  cases  (1)  the  volumetric  damage 
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dominates  the  deviatoric  one,  (2)  both  damage  modes  have  the  same  levels  of 
contributions  and  (3)  the  deviatoric  damage  mode  dominates  the  volumetric  one.  They 
are  compared  with  each  other. 

In  Figures  27  (a)  and  (b),  the  effective  stress-strain  curves  when  model  #1  and  #2  are 
used  are  presented,  respectively.  It  is  seen  that  the  macroscopic  behavior  of  both  the 
models  are  almost  the  same.  In  Figures  28-33,  the  distributions  of  stress  and  the  damage 
parameters  for  model  #1  are  depicted.  The  stress  distributions  at  10%  of  strain  are 
shown  in  Figure  28  (a)-(d).  The  distributions  look  alike.  However,  the  distributions  of 
the  volumetric  damage  parameters  are  quite  different  from  each  other,  as  shown  in 
Figures  29  (a)-(d).  The  value  of  damage  parameter  for  a  =  10.0  is  very  small.  For  the 
cases  of  a  =  0.1  and  a  =  0.99,  regions  with  relatively  large  volumetric  damage 
parameter  dv  run  perpendicular  to  the  loading  axis.  In  Figures  30  (a)-(d),  the 
distributions  of  the  deviatoric  damage  parameter  dD ,  when  the  block  is  subject  to  10% 
of  tensile  strain,  are  depicted.  The  damage  parameter  throughout  the  block  is  very  small 
as  shown  in  Figure  30  (b)  for  the  choice  of  parameter  a  =  0.1 .  The  locations/shapes  of 
relatively  low  and  high  damages  are  the  same  for  all  the  cases.  However,  it  is  clear  that, 
for  a  =  10.0 ,  the  regions  of  large  damage  parameter  run  in  45  degree  direction  form  the 
loading  axis.  The  trends  are  similar  to  the  case  of  isotropic  damage.  However,  the 
separate  volumetric/deviatoric  damage  model  with  a  =  10.0  gives  more  pronounced 
result. 

Figures  31  (a)-(d)  show  the  distributions  of  tensile  stress  when  the  block  is  subject  to 
15%  of  strain.  For  all  the  cases  in  Figure  31  (a)-(d),  the  regions  of  high  and  low  stresses 
run  in  the  vertical  direction.  The  concentrations  of  the  stress  are  more  pronounced  in  the 
case  of  a  =  10.0  than  the  other  cases.  The  distributions  in  Figures  31  (a),  (b)  and  (c) 
look  very  similar.  In  Figures  32  (a)-(d),  the  distributions  of  the  volumetric  damage 
parameter  are  depicted.  Volumetric  damage  regions  run  in  the  horizontal  direction  when 
a  =  0.1  is  assumed.  For  the  cases  of  a  =  0.99  and  a  =  10.0 ,  the  distributions  look  more 
uniform  than  that  of  a  =  0.1 .  There  are  localized  areas  with  relatively  small  value  of  the 
damage  parameter.  When  a  =  10.0  is  adopted,  the  value  of  the  damage  parameter  is 
small  throughout  the  block.  In  Figures  33  (a)-(d),  the  distributions  of  deviatoric  damage 
parameter  are  depicted.  As  expected,  the  value  is  small  throughout  the  block  when 
a  =  0.1 .  In  the  cases  of  a  =  0.99  and  a  =  10.0 ,  the  distributions  are  similar  to  each  other. 
There  are  the  bands  of  large  value  of  the  damage  parameter  that  run  along  45  degree 
direction  from  the  loading  axis.  The  directions  of  the  bands  are  slightly  different  from 
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the  case  of  the  isotropic  damage  model. 


From  the  results  presented  in  this  section,  it  is  clear  that  the  bands  of  large  volumetric 
damage  parameter  run  in  perpendicular  to  the  loading  direction,  whereas  those  of  large 
deviatoric  damage  parameter  are  along  the  45  degree  direction  from  the  tensile  axis.  The 
results  also  show  that  the  constant  a  can  control  which  damage  mode,  volumetric  or 
deviatoric,  dominates  the  other.  In  Figures  34-39,  the  results  obtained  based  on  the 
model  #2  are  shown.  The  overall  trends  are  the  same  as  those  obtained  based  on  the 
model  #1. 
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Figure  27  Effective  stress-strain  behaviors  that  were  generated  by  the  damage 
constitutive  models  with  the  material  stiffness  distribution  models  #1  and  #2. 
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Figure  28  Tensile  stress  distributions  at  10%  overall  tensile  strain  [(a)  isotropic 
damage  model,  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1, 
(c)  with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 
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Figure  29  The  distributions  of  dilatational  damage  parameter  at  10%  overall 
tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage  parameter),  (b) 
separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c)  with  a  =  0.99 
and  (d)  with  a  =  10.0  ]. 
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(a)  Isotropic  damage  ^  a  =  0  1 
model 


(c)  a  =  0.99 


(d)  a  =  10.0 


Figure  30  The  distributions  of  deviatoric  damage  parameter  at  10%  overall 
tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage  parameter),  (b) 
separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c)  with  a  =  0.99 
and  (d)  with  a  =  10.0  ]. 
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Figure  31  Tensile  stress  distributions  at  15%  overall  tensile  strain  [(a)  isotropic 
damage  model,  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1, 
(c)  with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 
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Figure  32  The  distributions  of  dilatational  damage  parameter  at  15%  overall 
tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage  parameter),  (b) 
separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c)  with  a  =  0.99 
and  (d)  with  a  =  10.0  ]. 


(a)  Isotropic  damage  (b)  a  =  0  1 
model 


(c)  a  =  0.99 


(d)  a  =  10.0 


Figure  33  The  distributions  of  deviatoric  damage  parameter  at  15%  overall 
tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage  parameter),  (b) 
separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c)  with  a  =  0.99 
and  (d)  with  a  =  10.0  ]. 
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Figure  34  Tensile  stress  distributions  at  10%  overall  tensile  strain  [(a)  isotropic 
damage  model,  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1, 
(c)  with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 
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Figure  35  The  distributions  of  dilatational  damage  parameter  at  10%  overall 
tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage  parameter),  (b) 
separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c)  with  a  =  0.99 
and  (d)  with  a  =  10.0  ]. 
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Figure  36  The  distributions  of  deviatoric  damage  parameter  at  10%  overall 
tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage  parameter),  (b) 
separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c)  with  a  =  0.99 
and  (d)  with  a  =  10.0  ]. 
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Figure  37  Tensile  stress  distributions  at  15%  overall  tensile  strain  [(a)  isotropic 
damage  model,  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1, 
(c)  with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 
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(a)  Isotropic  damage  (b)  a=0.l  (c )  a  =  0.99  (d)  a  =  10.0 

model 

Figure  38  The  distributions  of  dilatational  damage  parameter  at  15%  overall 
tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage  parameter),  (b) 
separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c)  with  a  =  0.99 
and  tdl  with  a  =  10.0.1 
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Figure  39  The  distributions  of  deviatoric  damage  parameter  at  15%  overall 
tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage  parameter),  (b) 
separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c)  with  a  =  0.99 
and  (d)  with  a  =  10.0  ]. 


4.3  The  behavior  of  a  cracked  specimen  -  isotropic  damage  model 

In  this  section,  we  present  results  for  cracked  specimen.  Finite  element  model  is  shown 
in  Figure  21.  Very  fine  meshes  are  placed  at  the  vicinity  of  the  crack  tip  and  relatively 
small  meshes  are  distributed  throughout  the  analysis  domain  unlike  the  cases  of 
ordinary  crack  analyses.  In  present  analysis,  we  assume  the  spatial  distribution  of 
material  stiffness  and  therefore  it  must  be  represented  by  the  fine  mesh  discretization. 
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The  distributions  of  material  stiffness  are  assumed  as  shown  in  Figure  20. 

The  distributions  of  tensile  stress  for  the  material  distribution  models  #1  and  #2  and  for 
homogeneous  material  case  are  shown  in  Figures  40  (a),  (b)  and  (c),  respectively,  when 
the  cracked  blocks  are  subject  to  0.3  %  of  overall  tensile  strain.  In  this  loading  step, 
entire  analysis  region  has  not  undergone  any  nonlinear  deformation.  For  the  cases,  in 
which  the  material  distributions  are  assumed,  we  can  see  that  there  are  the  bands  of 
relatively  high  stress  regions  running  in  vertical  direction.  In  Figures  41  (a),  (b)  and  (c) 
and  42  (a),  (b)  and  (c),  the  distributions  of  tensile  stress  are  presented  for  5%  and  7.5  % 
of  overall  strain,  respectively.  The  bands  of  relatively  high  stress  are  seen  in  Figures  41 
(a)  and  (b)  and  42  (a)  and  (b).  Other  than  these  points,  we  do  not  see  much  difference 
between  the  distributed  material  stiffness  models  and  the  homogeneous  model.  In 
Figure  43,  the  distributions  of  the  damage  parameter  are  presented  for  5%  of  overall 
strain.  It  is  seen  in  Figures  43  (a)  and  (b)  that  damage  zones  develop  at  the  crack  tip  and 
their  shapes  seem  to  be  influenced  by  the  material  stiffness  distributions.  When  the 
material  is  assumed  to  be  homogeneous  (Figure  43  (c)),  the  value  of  damage  parameter 
gradually  decreases  as  distance  from  the  crack  tip  increases.  In  the  cases  of  the 
distributed  material  stiffness  models,  the  distributions  of  the  damage  parameter  are  more 
complex.  We  can  see  that,  apart  from  the  crack  tip  damage  region,  there  are  some 
scattered  spots  with  relatively  large  value  of  the  damage  parameter.  The  shapes  of  crack 
tip  damage  regions  are  influenced  by  the  distributed  material  stiffness.  In  Figure  44,  the 
distributions  of  the  damage  parameter  are  shown  for  overall  strain  7.5%.  In  the  cases  of 
the  distributed  material  stiffness  models,  there  are  many  scattered  spots  having  large 
values  of  the  damage  parameter.  The  shapes  of  crack  tip  damage  zones  are  somewhat 
irregular. 
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Figure  40  Tensile  stress  distributions  at  overall  tensile  strain  0.3%  for  (a)  the 
distributed  material  stiffness  model  #1,  (b)  model  #2  and  (c)  homogeneous 
material,  with  the  isotropic  damage  constitutive  model. 
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(b)  Model  #2  (c)  Homogeneous 

Figure  41  Tensile  stress  distributions  at  overall  tensile  strain  5%  for  (a)  the 
distributed  material  stiffness  model  #1,  (b)  model  #2  and  (c)  homogeneous 
material,  with  the  isotropic  damage  constitutive  model. 


(a)  Model  #1 
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(a)  Model  #1 


(b)  Model  #2 


(c)  Homogeneous 


Figure  42  Tensile  stress  distributions  at  overall  tensile  strain  7.5%  for  (a)  the 
distributed  material  stiffness  model  #1,  (b)  model  #2  and  (c)  homogeneous 
material,  with  the  isotropic  damage  constitutive  model. 


(a)  Model  #1  (b)  Model  #2  (c)  Homogeneous 

Figure  43  The  distributions  of  the  damage  parameter  at  overall  tensile  strain  5% 
for  (a)  the  distributed  material  stiffness  model  #1,  (b)  model  #2  and  (c) 
homogeneous  material,  with  the  isotropic  damage  constitutive  model. 
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(a)  Model  #1  (b)  Model  #2  (c)  Homogeneous 

Figure  44  The  distributions  of  the  damage  parameter  at  overall  tensile  strain  7.5% 
for  (a)  the  distributed  material  stiffness  model  #1,  (b)  model  #2  and  (c) 
homogeneous  material,  with  the  isotropic  damage  constitutive  model. 


4.4  The  behavior  of  a  cracked  specimen  -  separate  isotropic/deviatoric  damage 
model 

In  this  section,  the  results  of  the  separate  isotropic/deviatoric  damage  model  for  the 
cracked  block  subject  to  tensile  deformation  are  described.  The  constant  a  which 
controls  the  damage  mode  (volumetric,  deviatoric  or  both  of  them)  is  set  to  be  0.1,  0.99 
or  10.0. 

In  Figure  45,  the  distributions  of  tensile  stress  are  presented  for  the  distributed  material 
stiffness  model  #1  when  the  cracked  block  is  subject  to  5%  of  tensile  deformation. 
Figures  45  (a),  (b),  (c)  and  (d)  show  the  results  of  the  isotropic  damage  model  and  the 
separate  isotropic/deviatoric  damage  model  with  a  =  0.1,  a  =  0.99  and  a  =  10.0, 
respectively.  For  all  the  material  models,  we  do  not  find  much  difference. 

In  Figure  46,  the  distributions  of  the  volumetric  damage  parameter  are  presented.  For  a 
comparison  purpose,  the  result  for  the  isotropic  damage  model  is  also  shown  in  Figure 
46  (a).  In  Figures  46  (b),  (c)  and  (d),  the  results  for  a  =  0.1 ,  a  =  0.99  and  a  =  10.0  are 
depicted,  respectively.  The  shapes  of  crack  tip  volumetric  damage  regions  are  almost  the 
same  in  Figures  46  (b),  (c)  and  (d).  Their  levels  are  different  as  the  constant  a  is  set  to 
be  0.1,  0.99  or  10.0.  In  Figure  47,  the  distributions  of  the  deviatoric  damage  parameter 
are  shown.  For  a  comparison  purpose,  the  result  of  the  isotropic  damage  model  is  also 
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depicted.  The  shapes  of  crack  tip  damage  zones  are  similar  to  each  other,  including  the 
case  of  the  isotropic  model.  However,  it  is  seen  that  as  the  value  of  the  constant  a 
increases,  the  damage  zone  tends  to  extend  to  upward  and  downward  directions.  When 
a  is  set  to  be  10.0,  the  shape  of  the  damage  zone  is  similar  to  that  of  plastic  zone  of 
elastic-plastic  analysis. 

In  Figure  48,  the  distributions  of  tensile  stress  are  presented  for  7%  of  overall  strain.  We 
do  not  find  much  difference  between  the  cases  of  the  isotropic  model  and  the  separate 
volumetric/deviatoric  damage  model  with  a  =  0.1,  a  =  0.99  and  a  =  10.0.  In  Figures 
49  (b),  (c)  and  (d),  the  distributions  of  the  volumetric  damage  parameters  are  presented. 
The  case  of  the  isotropic  damage  is  also  depicted  in  Figure  49  (a),  for  a  comparison 
purpose.  The  shapes  of  crack  tip  volumetric  damage  zones  are  almost  the  same. 
However,  their  levels  are  different  from  each  other,  as  the  values  of  the  constant  a  are 
set  to  be  0.1,  0.99  and  10.0.  In  Figures  50  (b),  (c)  and  (d),  the  distributions  of  the 
deviatoric  damage  parameter  are  depicted.  The  shapes  of  crack  tip  deviatoric  damage 
zones  in  the  cases  of  a  =  0.99  and  a  =  10.0  are  similar.  The  damage  zones  extend  in 
upward  and  downward  directions  from  the  crack  tip.  In  the  case  of  a  =  0.1 ,  the  damage 
zone  extends  in  the  forward  direction  from  the  crack  tip.  However,  the  value  of  the 
damage  parameter  is  very  small  compared  with  those  of  a  =  0.99  and  a  =  10.0 

In  Figures  51-56,  the  results  that  were  obtained  by  the  distributed  material  stiffness 
model  #2  are  presented.  We  find  the  same  observations  as  the  case  of  the  model  #1. 


(a)  Isotropic  (b)  a  =0.1  (c )  a  =  0.99  (d)  a  =  10.0 

Figure  45  Tensile  stress  distributions  of  crack  model  #1  at  5%  overall  tensile 
strain  [(a)  isotropic  damage  model,  (b)  separate  dilatational/deviatoric  damage 
model  with  a  =  0.1 ,  (c)  with  a  =  0.99  and  (d)  with  a  =  10.0], 
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(a)  Isotropic  (b)  or  =0.1  (c)  a  =  0.99  (d)  a  =  10.0 

Figure  46  The  distributions  of  the  volumetric  damage  parameter  of  crack  model 
#1  at  5%  overall  tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage 
parameter),  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c) 
with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 
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Figure  47  The  distributions  of  the  deviatoric  damage  parameter  of  crack  model 
#1  at  5%  overall  tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage 
parameter),  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (C) 
with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 
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Figure  48  Tensile  stress  distributions  of  crack  model  #1  at  7%  overall  tensile 
strain  [(a)  isotropic  damage  model,  (b)  separate  dilatational/deviatoric  damage 
model  with  a  =  0.1,  (c)  with  a  =  0.99  and  (d)  with  a  =  10.0], 
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(a)  Isotropic  (b)  a  =  0.1  (c)  a  =  0.99  (d)  a  =  10.0 

Figure  49  The  distributions  of  the  volumetric  damage  parameter  of  crack  model 
#1  at  7%  overall  tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage 
parameter),  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c) 
with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 
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Figure  50  The  distributions  of  the  deviatoric  damage  parameter  of  crack  model 
#1  at  7%  overall  tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage 
parameter),  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c) 
with  a  =  0.99  and  (d)  with  a  =  10.0], 


(a)  Isotropic  (b)  a  =0.1  (c )  a  =  0.99  (d)  a  =  10.0 

Figure  51  Tensile  stress  distributions  of  crack  model  #2  at  5%  overall  tensile 
strain  [(a)  isotropic  damage  model,  (b)  separate  dilatational/deviatoric  damage 
model  with  a  =  0.1 ,  (c)  with  a  =  0.99  and  (d)  with  a  =  10.0], 
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(a)  Isotropic  (b)  or  =0.1  (c)  a  =  0.99  (d)  a  =  10.0 

Figure  52  The  distributions  of  the  volumetric  damage  parameter  of  crack  model 
#2  at  5%  overall  tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage 
parameter),  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c) 
with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 


(a)  Isotropic  (b)  a  =  0.1  (c)  a  =  0.99  (d)  a  =  10.0 


Figure  53  The  distributions  of  the  deviatoric  damage  parameter  of  crack  model 
#2  at  5%  overall  tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage 
parameter),  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c) 
with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 
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Figure  54  Tensile  stress  distributions  of  crack  model  #2  at  7%  overall  tensile 
strain  [(a)  isotropic  damage  model,  (b)  separate  dilatational/deviatoric  damage 
model  with  a  =  0.1,  (c)  with  a  =  0.99  and  (d)  with  a  =  10.0], 
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(a)  Isotropic  (b)  a  =  0.1  (c)  a  =  0.99  (d)  a  =  10.0 

Figure  55  The  distributions  of  the  volumetric  damage  parameter  of  crack  model 
#2  at  7%  overall  tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage 
parameter),  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c) 
with  a  =  0.99  and  (d)  with  a  =  10.0  ]. 


(a)  Isotropic  (b)  a  =  0.1  (c)  a  =  0.99  (d)  a  =  10.0 

Figure  56  The  distributions  of  the  deviatoric  damage  parameter  of  crack  model 
#2  at  7%  overall  tensile  strain  [(a)  isotropic  damage  model  (the  isotropic  damage 
parameter),  (b)  separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  (c) 
with  a  =  0.99  and  (d)  with  a  =  10.0], 
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5.  Remarks 


5.1  Summary  of  present  investigation 

In  this  report,  the  results  of  analyses  in  which  material  is  assumed  to  be  polymeric 
particulate  composites  and  their  spatial  distributions  are  not  homogeneous  are  presented. 
The  uneven  distributions  of  the  particles  are  assumed  to  influence  the  stiffness  of  the 
material.  The  material  undergoes  damages  due  to  dewetting  between  the  particles  and 
the  polymeric  matrix  material  and  due  to  the  growth  of  micro  voids.  They  are  accounted 
for  by  the  damage  models.  In  order  to  model  the  uneven  distributions  of  the  particles, 
we  assumed  the  spatial  distributions  of  overall  material  stiffness. 

Two  types  of  damage  models  were  implemented  in  our  in-house  finite  element 
computer  program.  One  is  the  isotropic  damage  model  and  the  other  is  the  separate 
volume tric/deviatoric  damage  model.  Both  the  constitutive  models  are  based  on  the 
theory  of  Simo  and  Ju  [1],  The  separate  volumetric/de viatoric  model  was  developed 
during  this  course  of  study  in  order  to  model  material  damage  due  to  the  dilatational 
deformation  of  material.  This  damage  mode  is  mainly  due  to  the  growth  of  microvoids 
and  dewetting  between  the  particles  and  matrix  material. 

Since  the  material  properties  of  each  material  constituents  were  not  known,  we  fitted  the 
experimental  one-dimensional  stress-strain  curve  of  Kwon  and  Liu  [11]  by  the  damage 
constitutive  models.  Detailed  procedures  to  obtain  relationships  between  the  effective 
stress-like  parameters  and  the  damage  parameters  were  developed. 

The  damage  constitutive  models  (the  isotropic  and  the  separate  volumetric/de  viatoric 
damage  constitutive  models)  were  implemented  in  our  in-house  s-version  finite  element 
(s-FEM)  program.  The  s-FEM  computer  program  was  used  to  perform  the  analyses  of 
tensile  block  and  that  with  a  crack.  The  particles  in  the  polymeric  composites  distribute 
unevenly  at  a  local  level  but  evenly  at  larger  scale  level.  The  distribution  of  the  particles 
is  statistically  uniform  but  varies  locally.  Experimental  evidences  show  that 
representative  area  of  2  mm  x  2  mm  is  large  enough  so  that  the  particle  volume  fractions 
of  different  2  mm  x  2  mm  representative  areas  are  almost  the  same.  Thus,  we  built  two 
different  statistically  the  same  distributed  stiffness  models.  The  variations  of  the 
stiffness  of  the  material  reflect  the  variations  of  particle  volume  fraction. 
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Finite  element  analyses  were  performed  and  the  distributions  of  the  tensile  stress  and 
the  damage  parameters  were  visualized. 


5.2  Findings  from  the  analyses 

From  the  results  of  tension  blocks,  regions  of  high  and  low  stresses  connect  in  vertical 
direction.  The  distributions  at  the  same  overall  strain  are  similar  to  each  other  except  for 
the  case  of  a  =  10.0  at  15%  of  overall  strain.  The  differences  in  stress  values  at  high 
and  low  locations  are  more  pronounced  in  the  case  of  a  =  10.0  at  15%  of  overall  strain 
than  the  other  case.  It  is  considered  that  the  distributions  of  the  stress  are  governed  by 
the  spatial  variation  of  material  stiffness.  However,  the  values  are  influenced  by  how  the 
damage  is  sensitive  to  hydrostatic  or  deviatoric  stress.  That  is  why  a  =  10.0 ,  in  which 
the  material  damage  is  much  more  sensitive  to  the  deviatoric  stress  gives  somewhat 
different  result  from  the  others. 

The  distributions  of  damage  parameters  give  us  very  interesting  insights.  In  all  the  cases, 
we  find  the  bands  of  large  values  of  the  damage  parameters.  When  the  isotropic  damage 
model  is  adopted,  the  bands  tend  to  run  in  the  direction  40-50  degrees  from  the  loading 
axis.  When  the  separate  volumetric/deviatoric  damage  model  is  applied  to  the  analyses, 
the  bands  of  large  value  of  the  deviatoric  damage  parameter  run  in  about  45  degree 
direction  from  the  tensile  axis  and  those  of  the  volumetric  damage  parameter  are 
perpendicular  to  the  loading  direction.  This  means  that  the  bands  of  highly  damaged 
regions  are  influenced  not  only  by  the  spatial  variation  of  material  stiffness  but  also  by 
the  mode  of  damage.  In  the  case  of  a  =  10.0  in  which  the  deviatoric  damage  is 
considered  to  dominate  the  volumetric  mode,  the  direction  of  maximum  shear  stress. 
When  the  volumetric  damage  mode  dominates  the  other  ( a  =  0. 1 ),  the  damage  zones 
extend  in  the  section  of  maximum  tensile  stress.  It  is  considered  that  once  a  spot  of 
damaged  zone  is  present,  it  raises  the  tensile  stress  at  the  vicinity  of  the  weakened  spot 
in  the  direction  perpendicular  to  the  loading  axis,  as  shown  in  Figure  57.  Thus,  the 
hydrostatic  stress  at  that  position  enlarges,  resulting  in  the  progressive  material  damage. 
For  the  case  of  deviatoric  damage,  it  is  considered  that  the  damage  zones  tend  to  have 
the  45  degree  directions  from  the  loading  axis  for  the  similar  reason. 
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Figure  57  High  stress  spot  adjacent  to  the  damage  zone 


5.3  Relevance  to  the  case  of  polymeric  particulate  composite 

Though  in  this  research,  we  have  conducted  an  investigation  on  the  damage  behavior  of 
continuum  damage  solids.  The  growth  of  micro  voids  is  considered  to  be  governed  by 
the  hydrostatic  stress.  Therefore,  the  volumetric  damage  mode  takes  the  central  role  in 
such  a  case.  When  the  dewetting  between  the  particles  and  the  polymeric  matrix 
material  is  considered,  the  damage  mechanisms  may  be  more  complex  than  the  case  of 
the  growth  of  microvoids.  However,  after  a  particle  is  separated  from  matrix  material 
completely,  it  behaves  as  a  void.  Therefore,  the  growth  of  voids  plays  a  major  role  in  the 
deformation  mechanisms  of  the  composite. 

The  contributions  of  deviatoric  stress  components  to  the  material  damage  are  not  known. 
If  the  growth  of  micro  voids  and  the  dewetting  of  the  particles  were  the  major 
phenomena  in  the  damage  processes,  we  should  use  the  separate  volumetric/deviatoric 
damage  model  with  small  value  of  constant  a  . 
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6.  Conclusions 


In  this  investigation,  we  have  used  the  damage  constitutive  models  and  the  results  of 
different  models  have  been  compared  with  each  other.  Blocks  with  and  without  a  crack 
subject  to  tensile  deformation  were  analyzed  by  the  finite  element  method.  The  major 
findings  are: 

(1)  Bands  of  highly  damaged  zones  develop  in  the  tensile  block  and  their  directions 
depend  on  the  constitutive  model.  They  have  about  30,  45  and  90  degree  directions 
when  the  isotropic  and  the  separate  volume tric/deviatoric  damage  model  with  the 
constant  a  being  10  and  0.1,  respectively.  The  cases  of  a  =  10.0  and  a  =  0.1  almost 
correspond  to  those  of  deviatoric  and  dilatational  damage  models,  respectively. 

(2)  Two  different  statistically  the  same  material  stiffness  variations  in  the  tension  blocks 
were  assumed.  The  overall  trends  did  not  change  as  long  as  the  distributions  were 
statistically  the  same. 

(3)  The  shapes  of  crack  tip  damage  zones  depend  on  the  local  variations  of  the  material 
stiffness  at  the  crack  tip.  However,  the  overall  trends  were  the  same  for  the  same 
statistical  distributions  of  material  stiffness. 

(4)  The  shapes  of  the  crack  tip  damage  zone  depended  on  the  constitutive  models. 
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